home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / rng / ran3.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  2.9 KB  |  133 lines
  1. /* rng/ran3.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_rng.h>
  23.  
  24. /* This is an implementation of the algorithm used in Knuths's
  25.    subtractive generator, with the Numerical Recipe's ran3 paramters.
  26.    It is a subtractive lagged fibonnaci generator. */
  27.  
  28. static inline unsigned long int ran3_get (void *vstate);
  29. static double ran3_get_double (void *vstate);
  30. static void ran3_set (void *state, unsigned long int s);
  31.  
  32. #define M_BIG 1000000000
  33. #define M_SEED 161803398
  34.  
  35. typedef struct
  36.   {
  37.     unsigned int x;
  38.     unsigned int y;
  39.     unsigned long int buffer[56];
  40.   }
  41. ran3_state_t;
  42.  
  43. static inline unsigned long int
  44. ran3_get (void *vstate)
  45. {
  46.   ran3_state_t *state = (ran3_state_t *) vstate;
  47.   long int j;
  48.  
  49.   state->x++;
  50.  
  51.   if (state->x == 56)
  52.     state->x = 1;
  53.  
  54.   state->y++;
  55.  
  56.   if (state->y == 56)
  57.     state->y = 1;
  58.  
  59.   j = state->buffer[state->x] - state->buffer[state->y];
  60.  
  61.   if (j < 0)
  62.     j += M_BIG;
  63.  
  64.   state->buffer[state->x] = j;
  65.  
  66.   return j;
  67. }
  68.  
  69. static double
  70. ran3_get_double (void *vstate)
  71. {
  72.   return ran3_get (vstate) / (double) M_BIG ;
  73. }
  74.  
  75. static void
  76. ran3_set (void *vstate, unsigned long int s)
  77. {
  78.   ran3_state_t *state = (ran3_state_t *) vstate;
  79.   int i, i1;
  80.   long int j, k;
  81.  
  82.   if (s == 0)
  83.     s = 1;    /* default seed is 1 */
  84.  
  85.   j = (M_SEED - s) % M_BIG;
  86.  
  87.   /* the zeroth element is never used, but we initialize it for
  88.      consistency between states */
  89.  
  90.   state->buffer[0] = 0; 
  91.  
  92.   state->buffer[55] = j;
  93.  
  94.   k = 1;
  95.   for (i = 1; i < 55; i++)
  96.     {
  97.       int n = (21 * i) % 55;
  98.       state->buffer[n] = k;
  99.       k = j - k;
  100.       if (k < 0)
  101.     k += M_BIG;
  102.       j = state->buffer[n];
  103.  
  104.     }
  105.  
  106.   for (i1 = 0; i1 < 4; i1++)
  107.     {
  108.       for (i = 1; i < 56; i++)
  109.     {
  110.       long int t = state->buffer[i] - state->buffer[1 + (i + 30) % 55];
  111.       if (t < 0)
  112.         t += M_BIG;
  113.       state->buffer[i] = t;
  114.     }
  115.     }
  116.  
  117.   state->x = 0;
  118.   state->y = 31;
  119.  
  120.   return;
  121. }
  122.  
  123. static const gsl_rng_type ran3_type =
  124. {"ran3",            /* name */
  125.  M_BIG,                /* RAND_MAX */
  126.  0,                /* RAND_MIN */
  127.  sizeof (ran3_state_t),
  128.  &ran3_set,
  129.  &ran3_get,
  130.  &ran3_get_double};
  131.  
  132. const gsl_rng_type *gsl_rng_ran3 = &ran3_type;
  133.